import numpy as np
import scipy.io as scio
import copy


def rateAndDist(Ps, Zs, Hs):
    # print('p', Ps)
    m = len(Ps)
    n = len(Hs)
    Ps = np.array(Ps, dtype=int)
    Zs = np.array(Zs, dtype=int)
    C = np.zeros((1, n), dtype='double')
    PK = np.zeros((1, n), dtype='double')
    sig = np.sign(Zs - Ps)
    PK[0][Ps] = sig
    # print(PK,Ps,sig,Zs)
    for i in range(m):
        # print(Ps[i] + sig[i], Zs[i] - sig[i] + 1, sig[i])
        tem = np.arange(Ps[i] + sig[i], Zs[i], sig[i])
        C[0, tem] += sig[i]
    rates = sum(Hs[Ps])
    # print(C)
    # print(rates)
    Hs = Hs.reshape(511, -1)
    dists = 0.5 * (np.sum(np.dot((np.square(C) + np.square(C + PK)), Hs)))
    # print(rates, dists)
    return rates, dists


def caculateRD(Vs, maxLs, Gx, Hs):
    Rate = np.zeros(Gx)  # # 每个样本的容量与误差
    Dist = np.ones(Gx) * 100000000
    Layers = np.zeros(Gx)  # # 每个样本用到的峰值点和零值点对数
    hislength = np.shape(Hs)[0]
    for i in range(Gx):
        t_R = np.zeros(hislength)  # # 样本中每个直方图的容量与误差
        t_Ds = np.zeros(hislength)
        tmp_V = Vs[:, i * maxLs * 3: (i + 1) * maxLs * 3]
        t_Layer = np.zeros(hislength)  # # 样本中每个直方图用到的零点峰值点的对数
        for j in range(hislength):
            user = np.sum(tmp_V[j, maxLs * 2:maxLs * 3] == 1)
            if user == 0:
                t_R[j] = 0
                t_Ds[j] = 0
                continue
            User = tmp_V[j, maxLs * 2:maxLs * 3]
            Pset = tmp_V[j, 0:maxLs]
            Zset = tmp_V[j, maxLs:2 * maxLs]
            P = np.zeros(user)
            tmp_data = 0

            for z in range(maxLs):
                if User[z] == 1:
                    P[tmp_data] = Pset[z]
                    tmp_data = tmp_data + 1

            P_N_number = np.sum(P < 0)  # # 实际使用小于零点的峰值点个数
            P_P_number = np.sum(P >= 0)  # # 实际使用大于等于零点的峰值点个数
            Z_N_number = np.sum(Zset < 0)  # # 在零值点系列中小于零点的个数
            if P_N_number == 0:
                Z = Zset[Z_N_number:Z_N_number + P_P_number]
            elif P_P_number == 0:
                Z = Zset[Z_N_number - P_N_number:Z_N_number]
            else:
                NZ = Zset[Z_N_number - P_N_number:Z_N_number]
                PZ = Zset[Z_N_number:Z_N_number + P_P_number]
                Z = np.append(NZ, PZ)
            t_Layer[j] = len(P)
            P = P + 255
            Z = Z + 255
            t_R[j], t_Ds[j] = rateAndDist(P, Z, Hs[j, :])

        # if np.sum(t_Ds) != 0:
        Rate[i] = np.sum(t_R)
        Dist[i] = np.sum(t_Ds)
        Layers[i] = np.sum(t_Layer)
    return Rate, Dist, Layers


def find_best(R, D, G, rate, maxL, V):
    feature = np.zeros(G)
    for i in range(G):
        feature[i] = (R[i] >= rate)
    if np.sum(feature == 1) == 0:
        Map_PZ = []
        MP = []
        S_Rate = []
        S_dist = []
        success = 0
    t_D = D.copy()  ##寻找最好的样本
    ##容量不满足的情况下误差全部赋予inf
    for s in range(G):
        if feature[s] != 1:
            t_D[s] = 100000000
    pos = np.argmin(t_D)
    print('pos',pos)
    Vtest = V.copy()
    Best_V1 = Vtest[:, pos * 3 * maxL: (pos + 1) * 3 * maxL]
    Best_Rate1 = R[pos]
    Best_dist1 = D[pos]
    print(Best_dist1)
    ##对于误差的情况将样本分为两种情况，误差较大的feature标记为2，误差较小的标记为1.
    tmpdata = np.sum(t_D == 100000000)
    if G - tmpdata >= 20:  ##记录下容量满足的情况且误差在最大的前15个
        # pos = np.argsort(t_D, kind='mergesort')
        pos = np.argsort(t_D)
        BigD_pos = pos[G - tmpdata - 15:G - tmpdata]
        for j in range(len(BigD_pos)):
            feature[BigD_pos[j]] = 2
    return Best_V1, Best_dist1, feature


def promote_Rate(Vs, num, maxLs, Hs, USEs):
    ##容量不满足的情况下提升容量
    ##输入
    ##V整体样本
    ##num第几个样本
    ##峰值点零值点的最大点对数
    ##L直方图分配情况
    ##USE每个直方图候选的峰值点系列的上限
    ##输出
    ##V修改后的样本
    USEs = np.array(USEs, dtype=int)
    Vss = Vs.copy()
    tmp_V = Vss[:, num * maxLs * 3: (num + 1) * maxLs * 3]
    time = 1
    key = 0
    while not ((time > 3) or key):
        time = time + 1
        rand_num = np.random.randint(np.shape(Hs)[0])  ##随机选直方图的某一行
        # rand_num = 3
        use = tmp_V[rand_num, maxLs * 2:maxLs * 2 + USEs[rand_num]]
        p_set = tmp_V[rand_num, 0:USEs[rand_num]]  ##峰值点序列
        Nnumber = np.sum(p_set < 0)  ##峰值点序列中小于零的点的个数
        m = np.sum(use == 1)  ##找到所有使用峰值点的位置
        if len(use) != m:  ##在这一行中峰值点全都被使用，无法增加容量退出
            key = 1
    if key == 0:
        return Vs
    pose = np.where((use == 0) != 0)[0]  ##寻找所有没有使用的峰值点所在的位置
    # print('pose',pose)
    tmpdataa = len(pose)
    # print('tmpdata', tmpdataa)
    rand_num1 = np.random.randint(tmpdataa)
    # rand_num1 =0
    rand_num1 = pose[rand_num1]  ##随机寻找一个未使用的峰值点的位置
    # print('Nnumber', rand_num1, Nnumber)
    # print('use', use)

    if rand_num1 <= Nnumber - 1:  ##将选择的未使用的峰值点的位置向中心移动
        if rand_num1 == 0:  ##如果该点在零点左边的峰值点序列最左边直接使用它
            use[rand_num1] = 1
            tmp_V[rand_num, maxLs * 2: maxLs * 2 + USEs[rand_num]] = use

            Vs[:, num * maxLs * 3: (num + 1) * maxLs * 3] = tmp_V
            return Vs
        tmpdataa = np.sum(use[0:rand_num1] == 1)  ##随机选左边的点往右移
        if tmpdataa == 0:
            use[rand_num1] = 1
            tmp_V[rand_num, maxLs * 2: maxLs * 2 + USEs[rand_num]] = use
            Vs[:, num * maxLs * 3: (num + 1) * maxLs * 3] = tmp_V

        else:
            pose = np.where((use[0:rand_num1] == 1) != 0)[0]
            rand_num2 = np.random.randint(tmpdataa)
            rand_num2 = pose[rand_num2]
            use[rand_num2] = 0
            use[rand_num1] = 1

            tmp_V[rand_num, maxLs * 2: maxLs * 2 + USEs[rand_num]] = use
            Vs[:, num * maxLs * 3: (num + 1) * maxLs * 3] = tmp_V
    if rand_num1 > Nnumber - 1:
        if rand_num1 == USEs[rand_num] - 1:  ##如果该点在零点右边的峰值点序列最右边直接使用它
            use[rand_num1] = 1
            tmp_V[rand_num, maxLs * 2: maxLs * 2 + USEs[rand_num]] = use
            Vs[:, num * maxLs * 3: (num + 1) * maxLs * 3] = tmp_V
            return Vs
        tmpdataa = np.sum(use[rand_num1:USEs[rand_num]] == 1)  ##随机选右边的点往左移
        if tmpdataa == 0:
            use[rand_num1] = 1
            tmp_V[rand_num, maxLs * 2: maxLs * 2 + USEs[rand_num]] = use
            Vs[:, num * maxLs * 3: (num + 1) * maxLs * 3] = tmp_V
        else:
            pose = np.where((use[rand_num1:USEs[rand_num]] == 1) != 0)[0]
            rand_num2 = np.random.randint(tmpdataa)
            # rand_num2 =2
            rand_num2 = rand_num1 + pose[rand_num2]
            use[rand_num2] = 0
            use[rand_num1] = 1
            tmp_V[rand_num, maxLs * 2: maxLs * 2 + USEs[rand_num]] = use
            Vs[:, num * maxLs * 3: (num + 1) * maxLs * 3] = tmp_V
    return Vs


def descend_dist(Vs, num, maxLs, Hs, USEs):
    # # 容量满足的情况下减小误差
    # # 输入
    # # V整体样本
    # # num第几个样本
    # # 峰值点零值点的最大点对数
    # # L直方图分配情况
    # # USE每个直方图候选的峰值点系列的上限
    # # 输出
    # # V修改后的样本
    V = Vs.copy()
    tmp_V = V[:, num * maxLs * 3: (num + 1) * maxLs * 3]
    time = 1
    key = 0
    USEs = np.array(USEs, dtype=int)
    while not ((time > 3) or key):
        time = time + 1
        rand_num = int(np.random.randint(np.shape(Hs)[0]))  ##随机选直方图的某一行
        # rand_num = 4
        use = tmp_V[rand_num, maxLs * 2:maxLs * 2 + int(USEs[rand_num])]
        p_set = tmp_V[rand_num, 0:int(USEs[rand_num])]  ##峰值点序列
        Nnumber = np.sum(p_set < 0)  ##峰值点序列中小于零的点的个数
        m = np.sum(use == 0)  ##在这一行中峰值点没有被使用，无法衰减误差退出
        if len(use) != m:  ##在这一行中峰值点全都被使用，无法增加容量退出
            key = 1
    if key == 0:
        return Vs
    # print('use1', use)
    pose = np.where((use == 1) != 0)[0]  ##寻找所有使用的峰值点所在的位置
    tmpdatas = len(pose)
    # print('tmpdatas', tmpdatas)
    # print('pose', pose)
    rand_num1 = np.random.randint(tmpdatas)
    # rand_num1 = 0
    rand_num1 = pose[rand_num1]  ##随机寻找一个使用的峰值点的位置
    # print(10086)
    # print('Nnumber', Nnumber)
    # print('random', rand_num1)
    if rand_num1 <= Nnumber - 1:  ##将选择的未使用的峰值点的位置向中心移动
        if rand_num1 == 0:  ##如果该点在零点左边的峰值点序列最左边直接放弃它
            use[rand_num1] = 0
            tmp_V[rand_num, maxLs * 2: maxLs * 2 + USEs[rand_num]] = use
            Vs[:, num * maxLs * 3: (num + 1) * maxLs * 3] = tmp_V
            return Vs
        tmpdatas = np.sum(use[0:rand_num1] == 0)  ##随机选左边的点往左移
        if tmpdatas == 0:
            use[rand_num1] = 0
            tmp_V[rand_num, maxLs * 2: maxLs * 2 + USEs[rand_num]] = use
            Vs[:, num * maxLs * 3: (num + 1) * maxLs * 3] = tmp_V
        else:
            pose = np.where((use[0:rand_num1] == 0) != 0)[0]
            rand_num2 = np.random.randint(tmpdatas)
            rand_num2 = pose[rand_num2]
            use[rand_num2] = 1
            use[rand_num1] = 0
            tmp_V[rand_num, maxLs * 2: maxLs * 2 + USEs[rand_num]] = use
            Vs[:, num * maxLs * 3: (num + 1) * maxLs * 3] = tmp_V
    if rand_num1 > Nnumber - 1:
        if rand_num1 == USEs[rand_num] - 1:  ##如果该点在零点右边的峰值点序列最右边直接使用它
            use[rand_num1] = 0
            tmp_V[rand_num, maxLs * 2: maxLs * 2 + USEs[rand_num]] = use
            Vs[:, num * maxLs * 3: (num + 1) * maxLs * 3] = tmp_V
            return Vs
        tmpdatas = np.sum(use[rand_num1:int(USEs[rand_num])] == 0)  ##随机选右边的点往左移
        # print('2',tmpdatas)
        if tmpdatas == 0:
            use[rand_num1] = 0
            tmp_V[rand_num, maxLs * 2: maxLs * 2 + USEs[rand_num]] = use
            Vs[:, num * maxLs * 3: (num + 1) * maxLs * 3] = tmp_V
        else:
            pose = np.where((use[rand_num1:USEs[rand_num]] == 0) != 0)[0]
            # print('pse2',pose)
            rand_num2 = np.random.randint(tmpdatas)
            # rand_num2 =0
            rand_num2 = rand_num1 + pose[rand_num2]
            use[rand_num2] = 1
            use[rand_num1] = 0
            tmp_V[rand_num, maxLs * 2: maxLs * 2 + USEs[rand_num]] = use
            Vs[:, num * maxLs * 3: (num + 1) * maxLs * 3] = tmp_V
    return Vs


def descend_dist2(Vs, num, Best_V1s, maxLs, USEs):
    ##对于容量满足但是误差过大的情况直接使用历史最优样本（可能分割为16个直方图）中的十个子直方图的PZ点代替这些直方图中的相应位置
    # Vss = Vs.copy()
    tmp_V = Vs[:, num * maxLs * 3:(num + 1) * maxLs * 3]
    USEs = np.array(USEs, dtype=int)
    if np.shape(Vs)[0] < 10:  ##分割直方图的数量小于6个直接放弃
        mi = np.shape(Vs)[0]
    else:
        mi = 10
    rand_num = np.arange(np.shape(Vs)[0])
    np.random.shuffle(rand_num)  ##随机挑选十行

    rand_num = rand_num[0:mi]
    # print('rand_num', rand_num)
    # rand_num = np.array([2, 1, 10, 5, 6, 12, 8,11,9,15])
    for ms in range(len(rand_num)):
        use = Best_V1s[rand_num[ms], maxLs * 2: maxLs * 2 + USEs[rand_num[ms]]]
        # print('Best_V1', rand_num[k])
        # print('USE', USEs[rand_num[k]])
        # print('use.shape', use)
        # print('use',use.shape)
        # print('asdas', tmp_V[int(rand_num[k]), maxLs * 2: maxLs * 2 + USEs[rand_num[k]]].shape)
        tmp_V[int(rand_num[ms]), maxLs * 2: maxLs * 2 + USEs[rand_num[ms]]] = use
    Vs[:, num * maxLs * 3: (num + 1) * maxLs * 3] = tmp_V
    return Vs


def divide_rate(T, G, V, maxL, H, USE, rate, Best_V, feature, Best_dist):
    print('asd', Best_dist)
    final_bet = []
    for i in range(T):
        print('本轮best',Best_dist)
        for j in range(G):
            if feature[j] == 0:  # # 容量不满足的情况
                V = promote_Rate(V, j, maxL, H, USE)  # # 容量不满足的情况下提升容量
            elif feature[j] == 1:  # # 容量满足
                V = descend_dist(V, j, maxL, H, USE)  # # 容量满足的情况下减小误差
            elif feature[j] == 2:  # # 容量满足的情况下误差过大大规模减小误差
                V = descend_dist2(V, j, Best_V, maxL, USE)

        R, D, Layer = caculateRD(V, maxL, G, H)  # # 计算每个样本的容量与误差
        t_D = copy.deepcopy(D)
        feature = np.zeros(G)
        for z in range(G):
            # feature(z) = (R(z) >= search_rate + Layer(z) * 9)
            feature[z] = (R[z] >= rate)
        for k in range(G):
            if feature[k] != 1:
                t_D[k] = 1000000000
        # idx = np.where(t_D)
        pos = np.argmin(t_D)
        tmpdata = np.sum(t_D == 1000000000)
        if tmpdata != G:
            TMPDIST = D[pos]
            print('当前轮数',i)
            print('当前最佳',TMPDIST)
            print('历史最佳',Best_dist)
            if TMPDIST < Best_dist:
                Best_V = V[:, pos * 3 * maxL: (pos + 1) * 3 * maxL].copy()
                final_bet = Best_V.copy()
                # Best_Rate1 = R[pos]
                Best_dist = D[pos]
                print('当前最佳', TMPDIST)

        tmpdata = np.sum(t_D == 1000000000)
        if G - tmpdata > 20:  # # 记录下容量满足的情况且误差在最大的前15个
            # pos = np.argsort(t_D, kind='mergesort')
            pos = np.argsort(t_D)
            BigD_pos = pos[G - tmpdata - 15:G - tmpdata]
            for j in range(len(BigD_pos)):
                feature[BigD_pos[j]] = 2
    return final_bet, D, t_D, feature


def Choose_PZ(Best_V1, H, maxL):
    Map_PZ = np.zeros((np.shape(H)[0], maxL * 2 + 1))
    S_Rate = np.zeros(np.shape(H)[0])  ##样本中每个直方图的容量与误差
    S_dist = np.zeros(np.shape(H)[0])
    PzList = []

    for j in range(np.shape(H)[0]):
        user = np.sum(Best_V1[j, maxL * 2:maxL * 3] == 1)
        if user == 0:
            Map_PZ[j, 0:maxL * 2] = 255
            Map_PZ[j, -1] = 0
            PzList.append([[], []])
            continue
        User = Best_V1[j, maxL * 2:maxL * 3]
        Pset = Best_V1[j, 0:maxL]
        Zset = Best_V1[j, maxL:2 * maxL]
        P = np.zeros(user)
        tmp_data = 0
        for z in range(maxL):
            if User[z] == 1:
                P[tmp_data] = Pset[z]
                tmp_data = tmp_data + 1

        P_N_number = np.sum(P < 0)  ##实际使用小于零点的峰值点个数
        P_P_number = np.sum(P >= 0)  ##实际使用大于等于零点的峰值点个数
        Z_N_number = np.sum(Zset < 0)  ##在零值点系列中小于零点的个数
        if P_N_number == 0:
            Z = Zset[Z_N_number:Z_N_number + P_P_number]
        elif P_P_number == 0:
            Z = Zset[Z_N_number - P_N_number:Z_N_number]
        else:
            NZ = Zset[Z_N_number - P_N_number:Z_N_number]
            PZ = Zset[Z_N_number:Z_N_number + P_P_number]
            Z = np.append(NZ, PZ)
        layer = len(P)
        Map_PZ[j, 0:layer] = P
        Map_PZ[j, maxL:maxL + layer] = Z
        Map_PZ[j, -1] = layer
        P = P + 255
        Z = Z + 255
        PzList.append([P, Z])
        S_Rate[j], S_dist[j] = rateAndDist(P, Z, H[j, :])
    return S_Rate, S_dist, PzList


def initpopulation(Hs, maxL, G):
    # # V初始种群
    # # USE每个直方图实际上大于容量18的可用点的个数，如果可用点的对数大于maxL, 那么数值为maxL
    # # 就整个分割的直方图，构造初始值
    a = np.shape(Hs)[0]  # # 分割了多少个直方图
    Vs = np.zeros((a, maxL * 3 * G))
    # test = H[0, :].copy()
    # in_zero = np.where(test != 0)[0]

    USEs = np.zeros(a)  # # 每个直方图实际上大于容量18的可用点的个数，如果可用点的对数大于maxL, 那么数值为maxL
    for i in range(a):
        use_number = np.sum(Hs[i, :] > 100)
        if use_number > maxL:
            use_number = maxL

        pos = np.argsort(-Hs[i, :], kind='mergesort')
        # pos = np.argsort(-Hs[i, :])
        if i == 2:
            print('pos',pos)
            print('posnum',len(pos))
            print('use_number',use_number)
        # H[i,:] =-np.sort(-H[i,:])
        P = pos[0:use_number] - 255  # # 挑选峰值点
        P = np.sort(P)
        if i == 2:
            print('P',P)
        in_zero = np.where(Hs[i, :] != 0)[0]
        tmpdata = np.sum(P < 0)
        NZ = np.arange(in_zero[0] - 255 - tmpdata, in_zero[0] - 255)  # # 左边零值点
        tmpdata = np.sum(NZ < -255)
        if tmpdata > 1:
            print('crazyselect中左边零值点不够使用')

        tmpdata = np.sum(P >= 0)
        PZ = np.arange(in_zero[-1] - 255 + 1, in_zero[-1] - 255 + tmpdata + 1)  # # 右边非零点
        tmpdata = np.sum(PZ > 255)
        if tmpdata > 1:
            print('crazyselect中右边零值点不够使用')
        Z = np.append(NZ, PZ)  # # 挑选零值点
        layer = len(P)
        USEs[i] = layer
        Vs[i, 0: layer] = P
        if layer < maxL:
            Vs[i, layer: maxL] = 100000000

        Vs[i, maxL: maxL + layer] = Z
        if layer < maxL:
            Vs[i, maxL + layer: maxL * 2] = 100000000
        for j in range(1, G):
            Vs[i, j * maxL * 3: (j + 1) * maxL * 3] = Vs[i, 0: 3 * maxL]

    # V1 = V
    # save('V2.mat', 'V1')
    P_POS = np.zeros((a, 10))  # # 负的峰值点5个正的峰值点5个
    if a >= 8:
        temporary = 8
    else:
        temporary = a
    for i in range(temporary):
        P_set = Vs[i, 0:maxL]
        tmpdata1 = np.sum(P_set < 0)  # # 小于零的个数
        tmpdata2 = np.sum(P_set == 100000000)
        tmpdata2 = maxL - tmpdata1 - tmpdata2  # # 大于零的个数
        if tmpdata1 == 0:
            P_POS[i, 0: 5] = 100000000
        if tmpdata1 < 5:
            P_POS[i, 0: 5 - tmpdata1] = 100000000
            P_POS[i, 5 - tmpdata1: 5] = np.arange(1, tmpdata1 + 1)
        if tmpdata1 >= 5:
            P_POS[i, 0: 5] = np.arange(tmpdata1 - 4, tmpdata1 + 1)
        if tmpdata2 == 0:
            P_POS[i, 5: 10] = 100000000
        if tmpdata2 < 5:
            P_POS[i, 5 + tmpdata2: 10] = 100000000
            P_POS[i, 5: 5 + tmpdata2] = np.arange(tmpdata1 + 1, tmpdata1 + tmpdata2 + 1)
        if tmpdata2 >= 5:
            P_POS[i, 5: 10] = np.arange(tmpdata1 + 1, tmpdata1 + 5 + 1)
    if temporary == 8:
        P_POS[8: a, :] = 100000000
        tmpP_POS = np.array([
            P_POS[:, 5], P_POS[:, 4], P_POS[:, 6], P_POS[:, 3], P_POS[:, 7], P_POS[:, 2], P_POS[:, 8], P_POS[:, 1],
            P_POS[:, 9],
            P_POS[:, 0]]).T
        for i in range(10):  # # 对前面25个样本实行经验选值1
            tmpdata = tmpP_POS[:, 0: i + 1]
            tmpdata2 = np.shape(tmpdata)[1]
            for j in range(a):
                for z in range(tmpdata2):
                    if tmpdata[j, z] == 100000000:
                        continue
                    else:
                        # print('jz',tmpdata[j, z])
                        Vs[j, int(i * maxL * 3 + maxL * 2 + tmpdata[j, z] - 1)] = 1

    for i in range(10, 25):
        Vs[:, i * maxL * 3 + maxL * 2: i * maxL * 3 + maxL * 2 + maxL] = Vs[:,
                                                                         9 * maxL * 3 + maxL * 2:9 * maxL * 3 + maxL * 2 + maxL]

    for i in range(25, 28):  # # 对前面26至50个样本实行全选0
        for j in range(a):
            Vs[j, i * maxL * 3 + maxL * 2: i * maxL * 3 + maxL * 2 + int(USEs[j])] = 1

    for i in range(28, 50):  # # 对前面26至50个样本实行全选0
        for j in range(a):
            Vs[j, i * maxL * 3 + maxL * 2: i * maxL * 3 + maxL * 2 + int(USEs[j])] = 0
    # # 对51至100样本选取8对点分别1至4个直方图的0和 - 1
    # 点作为峰值点
    p11 = np.where((Vs[0, 0:maxL] == -1) != 0)[0]
    tmpdata = np.sum(Vs[0, 0:maxL] == -1)
    if tmpdata == 0:
        print('图像过于特殊')
    p12 = np.where((Vs[0, 0:maxL] == 0) != 0)[0]
    tmpdata = np.sum(Vs[0, 0:maxL] == 0)
    if tmpdata == 0:
        print('图像过于特殊')

    p21 = np.where((Vs[1, 0:maxL] == -1) != 0)[0]
    tmpdata = np.sum(Vs[1, 0:maxL] == -1)
    if tmpdata == 0:
        print('图像过于特殊')
    p22 = np.where((Vs[1, 0:maxL] == 0) != 0)[0]
    tmpdata = np.sum(Vs[1, 0:maxL] == 0)
    if tmpdata == 0:
        print('图像过于特殊')
    p31 = np.where((Vs[2, 0:maxL] == -1) != 0)[0]
    tmpdata = np.sum(Vs[2, 0:maxL] == -1)
    if tmpdata == 0:
        print('图像过于特殊')

    p32 = np.where((Vs[2, 0:maxL] == 0) != 0)[0]
    tmpdata = np.sum(Vs[2, 0:maxL] == 0)
    if tmpdata == 0:
        print('图像过于特殊')
    tmpdata = 50

    Vs[0, tmpdata * maxL * 3 + maxL * 2 + p11] = 1
    Vs[0, tmpdata * maxL * 3 + maxL * 2 + p12] = 1
    Vs[1, tmpdata * maxL * 3 + maxL * 2 + p21] = 1
    Vs[1, tmpdata * maxL * 3 + maxL * 2 + p22] = 1
    Vs[2, tmpdata * maxL * 3 + maxL * 2 + p31] = 1
    Vs[2, tmpdata * maxL * 3 + maxL * 2 + p32] = 1
    for i in range(51, 100):
        Vs[:, i * maxL * 3 + maxL * 2: (i + 1) * maxL * 3] = Vs[:,
                                                             tmpdata * maxL * 3 + maxL * 2: (tmpdata + 1) * maxL * 3]
    return Vs, USEs
